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Recent fluorescence spectroscopy measurements of single-enzyme kinetics have shown that enzy- 
matic turnovers form a renewal stochastic process in which the inverse of the mean waiting time 
between turnovers follows the Michaelis-Menten equation. Under typical physiological conditions, 
however, tens to thousands of enzymes react in catalyzing thousands to millions of substrates. We 
study enzyme kinetics at these physiologically relevant conditions through a master equation in- 
cluding stochasticity and molecular discreteness. From the exact solution of the master equation we 
find that the waiting times are neither independent nor are they identically distributed, implying 
that enzymatic turnovers form a non-renewal stochastic process. The inverse of the mean wait- 
ing time shows strong departures from the Michaelis-Menten equation. The waiting times between 
consecutive turnovers are anti-correlated, where short intervals are more likely to be followed by 
long intervals and vice versa. Correlations persist beyond consecutive turnovers indicating that 
multi-scale fluctuations govern enzyme kinetics. 

PACS numbers: 02.50.Ey,82.39.Fk,05.40.Ca 



Biological processes rely crucially on the catalytic ac- 
tivity of enzymes. In 1913, following the work of Wurtz 
and several others [1], Michaelis and Menten proposed [2] 
a reaction mechanism for catalysis where enzyme E binds 
reversibly with substrate S to form an enzyme-substrate 
complex ES which then dissociates irreversibly to form 

product P, while regenerating the enzyme: E + S s 1 ' 

k-i 

ES — ^> E + P. For thermodynamically large numbers 
of reactants, deterministic mass action kinetics provides 
the temporal variation of the concentrations of enzyme, 
complex and product. The rate of product formation 
is given by the classic Michaelis-Menten (MM) equation, 
provided suitable adiabaticity conditions are satisfied [3] . 

However, enzyme and substrate concentrations in bio- 
chemical catalysis are not thermodynamically large. In 
vivo enzyme concentrations vary from nanomolar to mi- 
cromolar, while the substrates are typically between a 
ten and ten thousand times more numerous [4]. An im- 
portant exception is in glycolysis where substrate con- 
centrations exceed those of enzymes [4] . In vitro enzyme 
concentrations vary from picomolar to nanomolar and 
substrates are typically a million times more numerous 
[5]. At these low concentrations, the inherent stochastic- 
ity of a single chemical reaction and the discrete change 
in the number of reactant molecules combine to gener- 
ate spontaneous, intrinsic fluctuations known as molecu- 
lar noise [6]. The temporal variation of catalysis, then, 
is also influenced by molecular noise and is a stochas- 
tic process in time. Recent advances in single molecule 
spectroscopy have been able to unravel some features of 
this stochastic process for catalysis involving a single en- 
zyme and numerous substrates [7, 8]. A striking feature 
is that the enzymatic turnovers generate a renewal point 



process where the waiting time t between product forma- 
tion events is independently and identically distributed. 
Remarkably, the inverse of the mean waiting time (t) -1 
obeys the MM equation which, in this interpretation, is 
valid not only for thermodynamically large systems, but 
also at the single-enzyme level. 

In this Letter, we study the stochastic process of 
enzymatic turnovers at concentrations between the ex- 
tremes of the thermodynamically large and single- 
enzyme regimes. In the thermodynamic limit the pro- 
cess reduces to deterministic evolution governed by mass 
action kinetics, while in the single-enzyme limit it re- 
duces to a renewal process. Our key findings are that 
for mesoscopic numbers of enzymes, the turnover pro- 
cess is of the non-renewal type with waiting times that 
are neither independent, nor identically distributed. We 
calculate the waiting time distributions and show that 
their inverse first moments do not obey the MM equa- 
tion. Consecutive waiting times are anti-correlated, with 
short intervals more likely to be followed by long inter- 
vals and vice-versa. The correlations persist beyond con- 
secutive turnovers and, depending on the number of en- 
zymes, can become substantially long-ranged. Together, 
these results imply that the enzymatic turnovers at the 
mesoscale cannot be described by mean production rates 
(as in the thermodynamic limit) or mean waiting times 
(as in the single-enzyme limit), but must be described 
by statistical measures which capture fluctuations over 
multiple time scales. 

Model. — We begin by describing the catalytic pro- 
cess through P(nEiTiESi n PTt)i the joint probability that 
there are tie enzymes, ues enzyme-substrate complexes 
and np products at any time t, starting initially with 
N enzymes, S substrates and no complexes or products. 
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Assuming that the system is well-mixed, the probability 
is taken to obey the Markovian chemical master equation 

P{n E ,n ESl n Pl t) = 
k a (n E + 1) P(n E + 1, n ES - 1, n P , t) 
+k_ 1 (n ES + l)P(n E - l,n ES + l,n P ,t) 
+k 2 {n ES + 1) P(n E - 1, n ES + 1, Hp - 1, t) 

k a n E + (fc-i + k 2 )n ES P(n E ,n ES ,n P ,t) (1) 

with the transition rates chosen to describe the MM 
kinetics of the enzyme catalysis reaction system given 
earlier. Since substrates are more numerous than en- 
zymes, the bimolecular second-order complexation step 

E + S ES is replaced by a pseudo-first order com- 
plexation step with an effective rate constant k a = k\S. 
The master equation generates stochastic trajectories of 
the kind shown in Fig. (1). Since enzymes are either 
converted to the enzyme-substrate complex or are re- 
generated from it, physical trajectories obey the con- 
straint n E + ties = N at all times. The probability 
distribution of these trajectories can then be written 
as P(N — n E s-,n E s,np,t\N, S) which we abbreviate to 
P(n E s,np, t\N). This simplifies the solution as there are 
two, and not three, independent variables. 

Exact solution. — We use the generating function 
method to obtain an exact solution of the master equa- 
tion. A related solution with n E and np as independent 
variables is given in [9] . Defining the generating function 
as 



G( Sl , S2 ,t) = j2J2 s " 



P(n ES ,n P ,t\N) (2) 



we find from the master equation its equation of motion 
[10], 

d t G(s u s 2 ,t) = k a N( Sl - l)G{ Sl ,s 2 ,t) + 

[(1 - Sl ){k b + k aSl ) - k 2 (l - 8i )] d sl G{ Sl ,s 2 ,t), (3) 

where k-\ + k 2 = k^. This partial differential equation in 
si,s 2 and t can be solved by the method of characteristics 
which, after a lengthy calculation, yields 



G(sx,s 2 ,t) 



1 



2 jv 

h - k a (l - 2si) 



-(B-A')t +e ~(B+A')t 



2A' 



-(B-A')t 



(B+A')t 



N 



(4) 



where A' = \y/{k a + h) 2 - 4fc a fc 2 (l - s 2 ) and B = 
\{k a + kb). P(n E s,np,t\N), obtained from the coeffi- 
cient of the Taylor expansion of the generating function 
in the si and s 2 variables, is an exact solution to the mas- 
ter equation. To support and complement this exact ana- 
lytical solution, we generate exact numerical trajectories 
of Eq. (1) using the Doob-Gillespie algorithm [11]. For 
the numerical simulations we non-dimcnsionalize time in 
units of k 2 and choose rate constants as k a = k 2 and 
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FIG. 1. (Color online) A trajectory of Eq. (1) for N = 2 
enzymes. The p-th product is generated at time T v . The 
waiting time between the p-th and (p + l)-th product is r v — 
2p+i — Tp. 



k-x = \k 2 . We generate ensembles of typically 10 6 tra- 
jectories to obtain the probability distributions of inter- 
est. One such trajectory is shown in Fig. (1). 

Turnover statistics. — The trajectories in Fig. (1) can 
be described in either of two ways : we can count the 
number of enzymatic turnover events np that have oc- 
curred in duration of time t, or, we can specify the time 
T p at which the p-th turnover occurs. The former, called 
the counting process description [12], has been used in 
most previous studies as it follows directly from the solu- 
tion P(n E s, np, t\N) of the master equation. The latter, 
called the point process description [13], has not (to the 
best of our knowledge) been studied before for multiple 
enzymes. This is the focus of our work. 

We define turnover times as T p = mi{t > : np(t) > 
p} for p = 1,2,..., which implies that T p < t if and 
only if np(t) > p. This provides the connection between 
the counting and point processes and relates the cumu- 
lative distribution of T p to that of n P by P(T p < t) = 
P(n p > p,t) [12]. Waiting times are defined from the 
turnover times by t p — T p — T p _i with the convention 
that To = 0. The point process is fully specified by the 
joint probability distributions w of cither the T p or the 
T p [12, 13]. Here we focus on the first-order distributions 
of the time to the p-th turnover w(T p ) and the interval 
between the p-th and (p+ l)-th turnovers w{t p ). We use 
second-order distributions w{T p ,T p+q ) to study correla- 
tions between the p-th and (p + ^)-th turnovers. 

First-order distributions. — We derive exact expres- 
sions for w(T p ) from the master equation solution us- 
ing the connection between the counting and point pro- 
cesses. It follows that P(T p < t) = P{n P > p,t) = 
l-P(n P <p,t) = 1-E^Lo^W)- Since P{n P ,t) = 
J2 n P{n E s,np,t\N), it follows by differentiation that 



>(Tp) = - E ^dtP(n E s,np,t\N) 



n„=0 n E s 



(5) 



t=T„ 
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FIG. 2. (Color online) Waiting time distributions for N — 
1000 and N = 1 enzymes. The waitings times are identi- 
cally distributed for a single enzyme (inset), but vary with 
the turnover number p for multiple enzymes. Solid lines are 
analytical results obtained from Eq. (7) while the symbols are 
simulation data. 
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FIG. 3. (Color online) First moment of T\ plotted, in 
Lineweaver-Burk fashion, against the inverse rate constant 
k a . Solid lines are first moments of Eq. (7) while symbols are 
simulation data. With more than one enzyme, the MM equa- 
tion (circles) is obeyed only in the limit of infinite substrate 
concentration or equivalently for l/k a — > 0. 



For fixed n p , the terms in the inner summation are related 
to the time derivative of G(s\, S2,t) evaluated at sx = 1 
through Eq. (2). The equation of motion, Eq. (3), is 
then used to eliminate the time derivative of G(s\, S2,t) 
in favour of its si derivative. From this, each term of the 
outer summation is obtained by taking rip derivatives 
with respect to s 2 , setting s 2 = and then summing 
over np to give [14] 

w(Tp)= (^l7! [arlasiG(si ' S2 ' Tp|iV)] -=M 2 =o- 

Using Eq. (4), w{T x ) = k 2 [d ai G(si,0,Ti\N)] ai=1 is ob- 
tained as 



k 2 k a N\ 



(2A) 



N 



e (A— B)Ti _ e -(A+B)T! 



(A + B)eS A - B ^ +(A- B)e-( B+A ^ 



N-l 



(7) 



and this is identical to w(ti). Here A = 
+ kb) 2 — 4fc a fc 2 and B = i(fc a + kb)- Since higher 
waiting times are differences of consecutive turnover 
times, the joint probability of consecutive turnovers is 
needed to calculate the w(r p ) for p > 1. This requires 
an involved calculation which we bypass by directly com- 
puting waiting time distributions from simulation trajec- 
tories. 

In Fig. (2) we compare w(t p ) for p = 1, 10, 100 for a 
reaction with N = 1000 enzymes. The t p are not identi- 
cally distributed. There is excellent agreement between 
the numerical result and the analytical calculation n, 
Eq. (7), plotted solid line. In the inset, we plot 
w(t p ), but for a single enzyme. The Tp are identically 
distributed, and agree with the analytical expression in 
Eq. (7) with N = 1. This clearly establishes the non- 
renewal nature of the turnover process when more than 
one enzyme participates in catalysis. 



Surprisingly, starting with a Markovian master equa- 
tion where waiting times between transitions are ex- 
ponentially distributed, we obtain a waiting time be- 
tween turnovers that is multi-exponential. For this, it 
is crucial to have more than one enzyme in the sys- 
tem. Then, as Eq. (5) shows, multiple internal states 
for the enzyme-substrate complex have to be summed 
over, and the resulting point process for the products 
alone is no longer Markovian. The multi-exponentiality 
of the waiting times is, therefore, consistent with the 
non-Markovian nature of the turnovers. For a single en- 
zyme, though, with only one internal enzyme-substrate 
state, there is no multi-exponentiality, but only a mono- 
exponential rise and fall. This is in agreement with ear- 
lier experimental [8, 15], numerical [7] and analytical [15] 
results. 

Moments of first- order distributions. — For a sin- 
gle enzyme, it follows from Eq. (7) that (Ti) = 
J °° dTx Tiw(Ti) = (S + K M )/k 2 S. The inverse of (T x ) 
then obeys the MM equation which has lead Xie and 
coworkers to extend the validity of the MM equation to 
the single-enzyme level [7, 8, 15]. However, for multiple 
enzymes, we find that the first moment no longer obeys 
the Michaelis-Menten equation as can be seen in in Fig. 
(3), where we plot N(T\) against l/k a in Lincweaver- 
Burk fashion. Thus, a turnover time interpretation of the 
MM equation is no longer valid for multiple enzymes. 

If N independent single-enzyme MM renewal process 
trajectories were to be pooled, there would be an iV-fold 
decrease in the mean turnover times. Fig. (3) shows that 
(7i) is larger than the MM estimate, indicating a slowing 
down of the kinetics due to cooperativity. The mean 
turnover time converges to the MM estimate only in the 
limit of infinite substrate concentration or equivalently 
for l/k a — > 0. Means (T p ) show similar behavior. 

The non-linearity in (Ti) is arises from the multi- 
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FIG. 4. The Pearson correlation coefficient of t p and fp+q as 
a function of lag q for TV — 10 (circles), TV = 50 (triangles) 
and TV = 100 (squares). The inset shows the joint distribution 
w(t p ,t p +i) for TV = 10. 



cxponentiality of w(Ti). For large N, there is no closed 
form analytical expression for the mean turnover time. 
However, in the limit of (k a + kf,) 2 ^> Ak a k 2 , which 
amounts to the steady-state approximation in the deter- 
ministic kinetics, the expression for the mean turnover 
time is given by 



1 



(TV - 1)TV 2 (5 6 



N 2 S 4 



k a k 2 (k a k 2 + (TV - 1)<5 2 ) 2 
(TV - 1)TV 2 5 6 



(k a k 2 + (TV - 1)<S 2 ) 2 ' k a k 2 (2k a k 2 + (TV - 2)<5 2 ) 2 J (8) 

This is obtained by using A — B w —5, A + B 2B and 
S = k a k 2 /(k a + kb) in Eq. (7), Taylor expanding to first 
order in (5/2B) exp [— (2B — S)Ti] and then computing 
the first moment of the approximated distribution [14]. 
The negative sign of the leading order correction term 
explains the curvature of the plot in Fig. (3). 

Second- order distributions and memory. — We com- 
pute the joint distributions w(t p , t p + q ) of t p and r p + q 
and their Pearson correlation coefficient from numerical 
trajectories. In Fig. (4) we plot the correlation coefficient 
against lag q, showing the joint distribution of consecu- 
tive intervals in the inset. The waiting times are anti- 
corrclated, where a short first interval is more likely to 
be followed by a long second interval and vice- versa. This 
memory effect shows a systematic variation with enzyme 
number, being strong and short-lived for fewer enzymes 
but weak and long-lived for more enzymes. With long- 
lived memory, fluctuation statistics will vary with the size 
of the temporal window, and multiple measures will be 
required to characterize the turnover process. In future 
work, we plan to explore this systematically, by studying 
higher-order joint distributions. The overall effect of the 
anti-correlations is to reduce the variance in the product 
turnovers when compared with a Poisson process. This 
may be biologically relevant to ensure a uniform rate of 
turnover in the steady state. 



Conclusion. — The non- renewal properties of enzymatic 
turnovers presented here can be verified by fluorescence 
experiments with well-mixed reactants. Fluctuations of 
intermediate states which must be summed over lead 
to multi-exponential waiting time distributions for the 
product and to the correlations between waiting times. 
These non-renewal aspects should appear in other mod- 
els of catalysis which involve several types of enzyme- 
substrate intermediates. Fluctuations of intermediate 
states can also provide a model for dynamic disorder, 
which has previously been modelled by fluctuating reac- 
tion rates. For second-order kinetics with substrate fluc- 
tuations [16], we numerically compute low order waiting 
time distributions and find negligible differences with our 
results. This justifies our use of pseudo-first order kinet- 
ics, which remains a reliable approximation at early times 
even when substrates fluctuate. In conclusion, the main 
implication of our work is that enzyme kinetics must be 
approached as a non-renewal stochastic process in time 
with fluctuations at multiple time scales. 
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